// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

/* 
 
 * NOTE: This file is the modified version of [s,d,c,z]panel_dfs.c file in SuperLU 
 
 * -- SuperLU routine (version 2.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * November 15, 1997
 *
 * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 *
 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 *
 * Permission is hereby granted to use or copy this program for any
 * purpose, provided the above notices are retained on all copies.
 * Permission to modify the code and to distribute modified code is
 * granted, provided the above notices are retained, and a notice that
 * the code was modified is included with the above copyright notice.
 */
#ifndef SPARSELU_PANEL_DFS_H
#define SPARSELU_PANEL_DFS_H

namespace Eigen {

namespace internal {

    template <typename IndexVector> struct panel_dfs_traits
    {
        typedef typename IndexVector::Scalar StorageIndex;
        panel_dfs_traits(Index jcol, StorageIndex* marker) : m_jcol(jcol), m_marker(marker) {}
        bool update_segrep(Index krep, StorageIndex jj)
        {
            if (m_marker[krep] < m_jcol)
            {
                m_marker[krep] = jj;
                return true;
            }
            return false;
        }
        void mem_expand(IndexVector& /*glu.lsub*/, Index /*nextl*/, Index /*chmark*/) {}
        enum
        {
            ExpandMem = false
        };
        Index m_jcol;
        StorageIndex* m_marker;
    };

    template <typename Scalar, typename StorageIndex>
    template <typename Traits>
    void SparseLUImpl<Scalar, StorageIndex>::dfs_kernel(const StorageIndex jj,
                                                        IndexVector& perm_r,
                                                        Index& nseg,
                                                        IndexVector& panel_lsub,
                                                        IndexVector& segrep,
                                                        Ref<IndexVector> repfnz_col,
                                                        IndexVector& xprune,
                                                        Ref<IndexVector> marker,
                                                        IndexVector& parent,
                                                        IndexVector& xplore,
                                                        GlobalLU_t& glu,
                                                        Index& nextl_col,
                                                        Index krow,
                                                        Traits& traits)
    {
        StorageIndex kmark = marker(krow);

        // For each unmarked krow of jj
        marker(krow) = jj;
        StorageIndex kperm = perm_r(krow);
        if (kperm == emptyIdxLU)
        {
            // krow is in L : place it in structure of L(*, jj)
            panel_lsub(nextl_col++) = StorageIndex(krow);  // krow is indexed into A

            traits.mem_expand(panel_lsub, nextl_col, kmark);
        }
        else
        {
            // krow is in U : if its supernode-representative krep
            // has been explored, update repfnz(*)
            // krep = supernode representative of the current row
            StorageIndex krep = glu.xsup(glu.supno(kperm) + 1) - 1;
            // First nonzero element in the current column:
            StorageIndex myfnz = repfnz_col(krep);

            if (myfnz != emptyIdxLU)
            {
                // Representative visited before
                if (myfnz > kperm)
                    repfnz_col(krep) = kperm;
            }
            else
            {
                // Otherwise, perform dfs starting at krep
                StorageIndex oldrep = emptyIdxLU;
                parent(krep) = oldrep;
                repfnz_col(krep) = kperm;
                StorageIndex xdfs = glu.xlsub(krep);
                Index maxdfs = xprune(krep);

                StorageIndex kpar;
                do
                {
                    // For each unmarked kchild of krep
                    while (xdfs < maxdfs)
                    {
                        StorageIndex kchild = glu.lsub(xdfs);
                        xdfs++;
                        StorageIndex chmark = marker(kchild);

                        if (chmark != jj)
                        {
                            marker(kchild) = jj;
                            StorageIndex chperm = perm_r(kchild);

                            if (chperm == emptyIdxLU)
                            {
                                // case kchild is in L: place it in L(*, j)
                                panel_lsub(nextl_col++) = kchild;
                                traits.mem_expand(panel_lsub, nextl_col, chmark);
                            }
                            else
                            {
                                // case kchild is in U :
                                // chrep = its supernode-rep. If its rep has been explored,
                                // update its repfnz(*)
                                StorageIndex chrep = glu.xsup(glu.supno(chperm) + 1) - 1;
                                myfnz = repfnz_col(chrep);

                                if (myfnz != emptyIdxLU)
                                {  // Visited before
                                    if (myfnz > chperm)
                                        repfnz_col(chrep) = chperm;
                                }
                                else
                                {  // Cont. dfs at snode-rep of kchild
                                    xplore(krep) = xdfs;
                                    oldrep = krep;
                                    krep = chrep;  // Go deeper down G(L)
                                    parent(krep) = oldrep;
                                    repfnz_col(krep) = chperm;
                                    xdfs = glu.xlsub(krep);
                                    maxdfs = xprune(krep);

                                }  // end if myfnz != -1
                            }      // end if chperm == -1

                        }  // end if chmark !=jj
                    }      // end while xdfs < maxdfs

                    // krow has no more unexplored nbrs :
                    //    Place snode-rep krep in postorder DFS, if this
                    //    segment is seen for the first time. (Note that
                    //    "repfnz(krep)" may change later.)
                    //    Baktrack dfs to its parent
                    if (traits.update_segrep(krep, jj))
                    //if (marker1(krep) < jcol )
                    {
                        segrep(nseg) = krep;
                        ++nseg;
                        //marker1(krep) = jj;
                    }

                    kpar = parent(krep);  // Pop recursion, mimic recursion
                    if (kpar == emptyIdxLU)
                        break;  // dfs done
                    krep = kpar;
                    xdfs = xplore(krep);
                    maxdfs = xprune(krep);

                } while (kpar != emptyIdxLU);  // Do until empty stack

            }  // end if (myfnz = -1)

        }  // end if (kperm == -1)
    }

    /**
 * \brief Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
 * 
 * A supernode representative is the last column of a supernode.
 * The nonzeros in U[*,j] are segments that end at supernodes representatives
 * 
 * The routine returns a list of the supernodal representatives 
 * in topological order of the dfs that generates them. This list is 
 * a superset of the topological order of each individual column within 
 * the panel.
 * The location of the first nonzero in each supernodal segment 
 * (supernodal entry location) is also returned. Each column has 
 * a separate list for this purpose. 
 * 
 * Two markers arrays are used for dfs :
 *    marker[i] == jj, if i was visited during dfs of current column jj;
 *    marker1[i] >= jcol, if i was visited by earlier columns in this panel; 
 * 
 * \param[in] m number of rows in the matrix
 * \param[in] w Panel size
 * \param[in] jcol Starting  column of the panel
 * \param[in] A Input matrix in column-major storage
 * \param[in] perm_r Row permutation
 * \param[out] nseg Number of U segments
 * \param[out] dense Accumulate the column vectors of the panel
 * \param[out] panel_lsub Subscripts of the row in the panel 
 * \param[out] segrep Segment representative i.e first nonzero row of each segment
 * \param[out] repfnz First nonzero location in each row
 * \param[out] xprune The pruned elimination tree
 * \param[out] marker work vector
 * \param  parent The elimination tree
 * \param xplore work vector
 * \param glu The global data structure
 * 
 */

    template <typename Scalar, typename StorageIndex>
    void SparseLUImpl<Scalar, StorageIndex>::panel_dfs(const Index m,
                                                       const Index w,
                                                       const Index jcol,
                                                       MatrixType& A,
                                                       IndexVector& perm_r,
                                                       Index& nseg,
                                                       ScalarVector& dense,
                                                       IndexVector& panel_lsub,
                                                       IndexVector& segrep,
                                                       IndexVector& repfnz,
                                                       IndexVector& xprune,
                                                       IndexVector& marker,
                                                       IndexVector& parent,
                                                       IndexVector& xplore,
                                                       GlobalLU_t& glu)
    {
        Index nextl_col;  // Next available position in panel_lsub[*,jj]

        // Initialize pointers
        VectorBlock<IndexVector> marker1(marker, m, m);
        nseg = 0;

        panel_dfs_traits<IndexVector> traits(jcol, marker1.data());

        // For each column in the panel
        for (StorageIndex jj = StorageIndex(jcol); jj < jcol + w; jj++)
        {
            nextl_col = (jj - jcol) * m;

            VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m);  // First nonzero location in each row
            VectorBlock<ScalarVector> dense_col(dense, nextl_col, m);   // Accumulate a column vector here

            // For each nnz in A[*, jj] do depth first search
            for (typename MatrixType::InnerIterator it(A, jj); it; ++it)
            {
                Index krow = it.row();
                dense_col(krow) = it.value();

                StorageIndex kmark = marker(krow);
                if (kmark == jj)
                    continue;  // krow visited before, go to the next nonzero

                dfs_kernel(jj, perm_r, nseg, panel_lsub, segrep, repfnz_col, xprune, marker, parent, xplore, glu, nextl_col, krow, traits);
            }  // end for nonzeros in column jj

        }  // end for column jj
    }

}  // end namespace internal
}  // end namespace Eigen

#endif  // SPARSELU_PANEL_DFS_H
